\documentclass[12pt, a4paper]{article}
\usepackage{geometry}
\usepackage{tikz}
\geometry{
	a4paper,
	left=12.7 mm,
	right=12.7 mm,
	top=12.7 mm,
	bottom=12.7 mm,
}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{physics}
\graphicspath{ {./} }
\usepackage{ctex}

\newcommand{\bvec}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\formula}[1]{\text{式} \ref{#1} }
\begin{document}
	\section{什么是微分方程，以及为什么要解微分方程}
	\footnote{参考：Sauer《数值分析》，欧阳洁等人的《数值分析》，MATLAB网站 https://ww2.mathworks.cn/help/matlab/math/choose-an-ode-solver.html}
	\textbf{微分方程}是含有函数导数项的方程。
	比如说，一个一阶常微分方程（ODE）形如
	\begin{equation}
		y = y(t) ~ \text{such that} ~ \dv{y}{t} = f(y,t)
	\end{equation}
	“一阶”指方程中函数导数的最高次为一；“常微分”指该函数只包括一个自变量。
	
	微分方程和代数方程有很大不同。
	求解微分方程$\dv{y}{t} = f(y,t)$，我们将得到一个\textbf{函数} $y = y(t)$；
	而求解代数方程（例如$x^2-2x+1=0$）我们会得到一个（或者一组）\textbf{数}。
	此外在求解微分方程时，只知道微分方程本身是不够的，
	我们还需要知道相应的初始条件等。
	
	在物理与工程中，不少动力学规律由微分方程描述。
	比如，我们熟悉的牛顿第二定律：
	$$
	m \dv[2]{\bvec r}{t} = \bvec F
	$$
	或者波动方程
	$$
	\pdv[2]{u}{t} - c^2 \pdv[2]{u}{x} = 0
	$$
	等都是微分方程。物理学的不少结论都蕴含在这些微分方程之中。
	因此，\textsl{学物理而不解微分方程，就如同喝可乐而不加冰、去旅游不发朋友圈}。
	
	\newpage
	\section{一阶ODE}
	\subsection{引入：放射性衰变}
	放射性元素的衰变方程是一个简单的一阶ODE：
	\begin{equation}
		N = N(t), 
		\left \{
		\begin{aligned}
			\dv{N(t)}{t} &= -\lambda N(t) \\
			N(t = t_0) &= N_0 \\
			t & \ge t_0 \\
		\end{aligned}
		\right.
	\end{equation}
	其中，$N$代表放射性原子的个数，$t$是时刻，而$\lambda$是一个常数。
	衰变方程的这种形式意味着，放射性元素的衰变速率和原子个数有关、原子越多衰变越快。
	
	接下来，我们想解出这个方程，以得到原子个数关于时间的关系$N=N(t)$。
	你可能已经知道了这个方程的解析解：
	$$
	N(t) = N_0 e^{-\lambda t} 
	$$
	但是对于大多数ODE方程，解析求解并不容易，并且经常是不可能的。
	因此，我们想直接数值求解它。
	所谓数值求解，就是求得一个个离散的时刻$t_0, t_1, t_2,..., $所对应的$N(t_0), N(t_1), N(t_2),...$。
	我们先假定，各时刻的间隔相同： $t_k = t_{k-1} +\Delta t = t_0 + k \Delta t, k=0,1,2,3,...$。
	
	\begin{table}[ht]
		\centering
		\caption{ODE数值求解}\label{tab_test251}
		\begin{tabular}{cccccc}
			\hline
			$t$ & $t_0$  & $t_1 = t_0+\Delta t$ & $t_2= t_1+\Delta t$ & ...\\
			\hline
			$N$ & $N(t_0) = N_0$ & $N(t_1)$ & $N(t_2)$ & ...\\
			\hline
		\end{tabular}
	\end{table}
	
	\newpage
	
	\subsection{显式Euler法}
	\begin{figure}[h]
		\centering
		\includegraphics[width=0.5\linewidth]{expiter}
		\caption{显Euler迭代示意图}
		\label{fig:expiter}
	\end{figure}
	
	显式Euler法可能是最简单的ODE数值解法之一，其直接源于对导数含义的原始理解：
	\begin{equation}
		\dv{N(t)}{t} = -\lambda N(t)  \Rightarrow \frac{N(t+\Delta t) - N(t)}{\Delta t} = -\lambda N(t) \Rightarrow N(t+\Delta t)  = N(t) -\lambda N(t) \Delta t 
	\end{equation}
	下一时刻的$N(t+\Delta t) $可由上一时刻$N(t) $的值推导得到。很少数学公式能有这么直观的形式！
	很容易改写为方便计算机实现的递推公式：
	\begin{equation}
		N^{(k+1)} = N^{(k)} - \lambda N^{(k)}\Delta t \qquad k = 0,1,2,3,...
	\end{equation}
	或（仅仅是书写的不同）
	\begin{equation}
		N^{(k)} = N^{(k-1)} - \lambda N^{(k-1)}\Delta t \qquad k = 1,2,3,...
	\end{equation}
	此处以$k$上标代表时刻以简化书写，$N^{(k)} = N(t_k) = N(t_0 + k \Delta t)$。
	
	对于一个一般性的ODE
	\begin{equation}
		y = y(t), 
		\left \{
		\begin{aligned}
			\dv{y(t)}{t} &= f(y, t) \\
			y(t_0) &= y_0 \\
			t & \ge t_0 \\
		\end{aligned}
		\right.
	\end{equation}
	其中$f$是已知函数。其显式Euler格式是
	\begin{equation}
		\boxed
		{
			y^{(k+1)} = y^{(k)} + f(y^{(k)}, t_k) \Delta t \qquad k = 0,1,2,...
		}
	\end{equation}
	
	显式Euler法非常直观清晰简洁明了，但有一个缺点：精度低、速度慢。
	如果你给的$\Delta t$太大、那迭代结果很快与实际的结果误差太大；
	为此你必须将$\Delta t$取得很小，而这又导致需要非常庞大的计算量！
	我们继续介绍一些改进方法、这些方法具有更高的精度，以及更复杂的数学形式。
	
	\newpage
	\subsection{Euler梯形法(Euler预估-校正公式)}
	
	\begin{figure}[h]
		\centering
		\includegraphics[width=0.5\linewidth]{trapiter}
		\caption{Euler预估-校正公式迭代示意图.橙色线与虚线圆代表预估过程}
		\label{fig:trapiter}
	\end{figure}
	
	Euler梯形法有助于提高精度。字面上看，Euler梯形法平均了上下两步的导数：
	$$
	y^{(k+1)} = y^{(k)} + \frac{f(y^{(k+1)}, t_{k+1}) + f(y^{(k)}, t_{k})}{2} \Delta t\\
	$$
	但细心的你很快发现，右侧$f(y^{(k+1)}, t_{i+1})$项中已经使用了待求的$y^{(k+1)}$。于是我们遇到了一个\textsl{先有蛋还是先有鸡的哲学问题}：如果我们还不知道$y^{(k+1)}$，那么我们怎么使用它来求解$y^{(k+1)}$？
	
	一种方法是，先用显式Euler法 “预估”$y^{(k+1)}$，再将其代入公式“校正”。总之：
	\begin{equation}
		\boxed
		{
			\begin{aligned}
				\overline{y^{(k+1)}} &= y^{(k)} + f(y^{(k)}, t_k) \Delta t  \\
				y^{(k+1)} &= y^{(k)} + \frac{f(\overline{y^{(k+1)}}, t_{k+1}) + f(y^{(k)}, t_{k})}{2} \Delta t\\
			\end{aligned}
			\qquad k =0,1,2,...
		}
	\end{equation}
	$\overline{y^{(k+1)}} $是“预估”。
	
	因为Euler法平均了两步的导数，有点像数值积分中的梯形公式，因此称为“Euler梯形公式”；
	又因为其用到了显式Euler的预估结果来校正结果，因此也称“Euler预估-校正公式”。
	
	\subsection{Taylor法}
	如果我们直接Taylor展开$y(t_k+\Delta t)$，我们会得到：
	\begin{equation}
		y(t_k+\Delta t) = y(t_k) + y'(t_k) \Delta t + \frac{1}{2} y''(t_k) (\Delta t)^2 + ...
	\end{equation}
	如果我们只取一阶项，
	$$
	y(t_k+\Delta t) = y(t_k) + y'(t_k) \Delta t 
	$$
	其实就是显式Euler法，其中$y'(t_k) = f(y^{(k)}, t_k)$。
	
	Taylor法的要义在于，我们还要取高阶项以获得更高精度。但高阶项，比如$y''(t_k)$，是什么？
	这就有点复杂了，事实上
	\begin{equation}
		y''(t_k) = \dv{y'(t_k)}{t} =  \dv{f(y^{(k)},t_k)}{t} = \pdv{f}{t}|_{y=y^{(k)}} \dv{y}{t}|_{t=t_k} + \pdv{f}{t}_{y=y^{(k)}} = f_y f +f_t
	\end{equation}
	其中$f_y$是$f$关于$y$的偏导、$f_t$是$f$关于$t$的偏导。
	高阶导数需要链式法则，并且需要知道$f$关于$y,t$的偏导数。
	随着阶数升高，这个公式会变得又臭又长。
	
	一个二阶的Taylor法是
	\begin{equation}
		\boxed
		{
			y^{(k+1)} = y^{(k)} + f(y^{(k)}, t_k) \Delta t + \frac{1}{2} \left(f_y(y^{(k)}, t_k) f(y^{(k)}, t_k) + f_t(y^{(k)}, t_k)\right) (\Delta t)^2 + ...
			\qquad k =0,1,2,...
		}
	\end{equation}
	
	\subsection{龙格-库塔法}
	龙格-库塔法是Taylor法的一种实现，但巧妙地避免了求偏导数。
	二阶龙格-库塔法是
	\begin{equation}
		\boxed
		{
			y^{(k+1)} = y^{(k)} 
			+ f\left(y^{(k)} + \frac{f(y^{(k)}, t_k) \Delta t}{2}, t_k + \frac{\Delta t}{2}\right)\Delta t
			\qquad k =0,1,2,...
		}
	\end{equation}
	为什么二阶龙格-库塔法是这样的呢？我们对这个公式做Taylor展开：
	$$
	\begin{aligned}
		y^{(k+1)} &= y^{(k)} 
		+ f\left(y^{(k)} + \frac{f(y^{(k)}, t_k) \Delta t}{2}, t_k + \frac{\Delta t}{2}\right)\Delta t \\
		& \approx y^{(k)} + \left( f(y^{(k)}, t_k) + \pdv{f}{y} \frac{f(y^{(k)}, t_k) \Delta t}{2} + \pdv{f}{t}  \frac{\Delta t}{2} \right) \Delta t\\
		& = y^{(k)} + f(y^{(k)}, t_k) \Delta t + \frac{1}{2}(f_y f + f_t)(\Delta t)^2
	\end{aligned}
	$$
	整理后发现其与二阶Taylor法是一样的。
	一般而言，四阶龙格-库塔法（术语：RK4）比较常用。
	数值分析教材一般会直接给出相应的公式，不需要\textsl{重复造轮子}。
	
	
	
	\newpage
	\section{ODE方程组}
	\subsection{引入：化学反应}
	假设我们有一个化学反应
	\begin{equation}
		A + B \to 2C
	\end{equation}
	当化学反应按这个式子发生一点时（化学反应进度变化$\dd \xi$后），
	由于物质守恒，各种物质的量将成比例变化，比例为化学计量数。
	化学计量数是化学物质前的系数，记生成物为正、反应物为负。
	在这个例子中，
	$A$的化学计量数 $\nu_A = -1$, 
	$B$的化学计量数 $\nu_B = -1$, 
	$C$的化学计量数 $\nu_C = 2$。
	\begin{equation}
		\text{发生一点化学反应 $\dd \xi$后：} 
		\qquad
		\dd c_A = \nu_A \dd \xi =  - \dd \xi, ~
		\dd c_B =  - \dd \xi, ~
		\dd c_C = 2 \dd \xi
	\end{equation}
	也就是说，（只考虑正向反应，）若化学反应速率为$v = \dv{\xi}{t}$，
	那么
	\begin{equation}
		\dv{c_A}{t} = -v, \dv{c_B}{t} = -v, \dv{c_C}{t} = 2v 
	\end{equation}
	此外，化学反应速率 $v$ 往往与各物质含量有关：
	\begin{equation}
		v = v (c_A, c_B, c_C, ...)
	\end{equation}
	因此，各物质含量关于时间的变化可以被表达为一组ODE方程组
	\begin{equation}
		c_A = c_A(t), c_B = c_B (t),  c_C = c_C (t)
		\qquad
		\left \{
		\begin{aligned}
			\dv{c_A}{t} & = -v (c_A, c_B, c_C, ...) \\
			\dv{c_B}{t} & = -v (c_A, c_B, c_C, ...) \\
			\dv{c_C}{t} & = 2v (c_A, c_B, c_C, ...) \\
			c_A(t_0) &= c_{A0} \\
			c_B(t_0) &= c_{B0} \\
			c_C(t_0) &= c_{C0} \\
		\end{aligned}
		\right.
	\end{equation}
	
	\subsection{显式Euler法}
	有时，我们会遇到一组互相关联的ODE方程：
	\begin{equation}
		y_1 = y_1(t), y_2 = y_2 (t), ...
		\qquad
		\left \{
		\begin{aligned}
			\dv{y_1}{t} & = f(y_1, y_2,...,t) \\
			\dv{y_2}{t} & = f(y_1, y_2,...,t) \\
			& ... \\
			y_1(t_0) & = (y_1)_0 \\
			y_2(t_0) & = (y_2)_0 \\
			& ... \\
		\end{aligned}
		\right.
	\end{equation}
	解决这个问题的方法也很简单：只需要分别对每一个$y_i$做显式Euler迭代：
	\begin{equation}
		\begin{aligned}
			y_1^{(k+1)} &= y_1^{(k)} + f(y_1^{(k)}, y_2^{(k)}, ..., t_k) \Delta t \\
			y_2^{(k+1)} &= y_2^{(k)} + f(y_1^{(k)}, y_2^{(k)}, ..., t_k) \Delta t \\
			& ...
		\end{aligned} 
		\qquad k=0,1,2,...
	\end{equation}
	这样，我们就能解出ODE方程组。上述表达式还可以压缩为向量形式：
	\begin{equation}
		\bvec y = (y_1,y_2,...)^T \qquad \bvec y^{(k+1)} = \bvec y^{(k)} + f(\bvec y^{(k)}, t_k) \Delta t
	\end{equation}
	形式上重新将其变为了熟悉的模样，同理可运用上述的龙格-库塔法等
	
	\section{高阶ODE}
	我们还会遇到高阶ODE，例如包含二阶导数的ODE方程：
	\begin{equation}
		y = y(t) \qquad 
		\left \{
		\begin{aligned}
			\dv[2]{y}{t} &= f(y,t) \\
			y(t_0) &= y_0 \\
			\dv{y}{t} |_{t=t_0} & = y'_0 \\
		\end{aligned}
		\right.
	\end{equation}
	解决高阶ODE的一种方法是，通过代换将其转换为低阶ODE方程组。设
	\begin{equation}
		v = \dv{y}{t}
	\end{equation}
	那么，二阶ODE降阶为两个一阶ODE：
	\begin{equation}
		\dv[2]{y}{t} = f(y,t)
		\Rightarrow 
		\left \{
		\begin{aligned}
			\dv{y}{t} & = v \\
			\dv{v}{t} & = f(y,t) \\
		\end{aligned}
		\right.
	\end{equation}
	这是我们已经会解的问题：
	\begin{equation} \label{eq_highorder_ode}
		\begin{aligned}
			y^{(k+1)} &= y^{(k)} + v^{(k)} \Delta t \\
			v^{(k+1)} & = v^{(k)} + f(y^{(k)}, t_k) \Delta t \\
		\end{aligned}
		\qquad k=0,1,2,...
	\end{equation}
	可见，为了迭代$y$与$v$，我们需要两个初始条件，$y^{(0)}$与$v^{(0)} = y'_0$。
	这就是为什么对于高阶ODE，我们不仅得知道函数的初始值，还得知道函数导数的初始值等。
	
	我们很想用这种方法求解牛顿第二定律----
	但是先别急，牛顿第二定律有点特殊，最好\textbf{不要}用这种方法求解。
	请参考隔壁Verlet算法。

\end{document}

